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ABSTRACT 

Highly  accurate  far  field  computational  boundary  conditions  for  inviscid^ 
two-dimensional  isentropic  duct  flow  problems  are  developed  from  analytic 
solutions  of  the  linearized,  second-order  Euler  equations.  The  Euler 
equations  are  linearized  about  a  constant  pressure,  rectilinear  flow 
condition.   The  boundary  procedure  can  be  used  with  any  numerical  Euler 
solution  method  and  allows  computational  boundaries  to  be  located  extremely 
close  to  the  nonlinear  region  of  interest.  Numerical  results  are  presented 
which  show  that  the  boundary  conditions  and  far  field  analytic  solutions 
provide  a  smooth  transition  across  a  computational  boundary  to  the  true  far 
field  conditions  at  infinity.  The  cost  of  upgrading  first-order  boundary 
conditions  to  second-order  is  slight. 


I.   INTRODUCTION 


Numerical  solution  procedures  for  nonlinear  fluid  dynamic  equations 
usually  use  one  or  more  artificial  computational  boundaries  located  at  some 
distance  from  the  primary  region  of  interest  in  order  to  limit  the  physical 
domain  to  finite  size.  If  the  flow  crossing  such  a  boundary  (either  inflow  or 
outflow)  is  subsonic,  then  some  type  of  computational  boundary  conditions  must 
be  imposed  which  simulate  the  influence  of  the  true  far  field  conditions  at 
infinity.  These  boundary  conditions  must  be  such  that  waves  crossing  the 
boundary  do  not  produce  erroneous  reflections  back  into  the  computational 
field  to  degrade  the  calculations.  It  is  generally  acknowledged  that  simply 
imposing  free  stream  conditions  (or  conditions  at  infinity)  at  computational 
boundaries  is  usually  inappropriate  because  of  the  spurious  reflections  back 
into  the  computational  domain  which  result.  Standard  practice  has  consisted 
of  locating  the  boundaries  quite  far  from  the  region  of  interest  in  an  attempt 
to  simplify  the  boundary  condition  models  and  minimize  any  effects  of  incon- 
sistent modeling.  The  net  effect  is  a  significant  increase  in  the  number  of 
grid  points  required  for  an  accurate  flowfield  calculation. 

A  boundary  modeling  procedure  for  two-dimensional  internal  flows  was 
presented  in  Reference  1  which  alleviates  the  difficulties  mentioned  above  and 
also  allows  the  computational  boundaries  to  be  located  much  closer  to  the 
nonlinear  region  of  interest.  The  procedure  is  limited  to  steady,  inviscid 
flow,  although  the  flow  can  be  rotational.   It  represents  a  logical  first- 
order  extension  of  the  so-called  characteristic  (or  zero-order)  boundary 
conditions  commonly  used  with  inviscid  numerical  solution  methods.  Extension 
to  axisymmetric  or  three-dimensional  flows  is  straightforward. 

The  analysis  presented  here  extends  the  first-order  analysis  to  second 
order  for  isentropic  duct  flow  and  provides  a  logical  extension,  in  an 
asymptotic  sense,  of  the  first-order  analysis.   It  also  illustrates  a  consis- 
tent procedure  for  coupling  linearized  analytic  solutions  with  nonlinear 
numerical  solutions  by  means  of  computational  boundary  conditions.  The 
greater  accuracy  of  the  second-order  boundary  procedure  allows  the  computa- 
tional boundaries  to  be  placed  even  closer  to  the  nonlinear  computational 


region  thereby  further  reducing  the  number  of  grid  points  needed  for  the 
numerical  solution.  Even  though  the  analysis  and  boundary  conditions  derived 
from  it  are  limited  to  isentropic  flow,  many  non-isentropic  duct  flow  problems 
have  isentropic  flow  crossing  the  inflow  boundary  at  which  these  second-order 
boundary  conditions  can  be  imposed. 

The  present  analysis  is  based  on  the  Riemann  variable  formulation  of  the 
Euler  equations  given  in  Reference  2   simplified  by  the  isentropic  assumption. 
This  represents  a  natural  starting  point  because  the  characteristic  (or 
zero-order)  boundary  conditions  mentioned  above  are  expressed  in  terms  of 
Riemann  variables.  The  equations  are  linearized  about  a  constant  pressure, 
rectilinear  flow  condition,  which  truly  represents  conditions  at  infinity. 
These  linearized  equations  are  assumed  applicable  in  the  far  field  region 
beyond  a  computational  boundary.  Within  the  nonlinear  computational  domain, 
strong  entropy-producing  (i.e.,  rotational)  effects  can  exist  which  create 
variations  in  density,  velocity,  etc.  in  the  downstream  far  field  in  the 
streamline-normal  direction  which  are  not  necessarily  small  perturbations. 
Such  variations  were  modeled  in  the  analysis  presented  in  Reference  1. 
However,  because  of  the  two-level  perturbation  procedure  and  subsequent 
approximate  solution  method  used  in  that  analysis,  inclusion  of  entropy 
effects  in  the  second-order  analysis  does  not  seem  justified.  As  mentioned 
above,  however,  the  second-order  isentropic  analysis  is  still  applicable  to 
non-isentropic  problems  having  isentropic  inflow  conditions. 

The  first-order  linearized  equations  analyzed  in  Reference  1  are  homo- 
geneous and  were  solved  using  separation  of  variables  and  Fourier  analysis. 
The  second-order  equations  treated  herein  are  non-homogeneous,  but  can  also  be 
solved  by  the  same  techniques.  These  equations  are  written  in  different  forms 
for  the  upstream  and  downstream  far  field  regions  because  the  downstream  flow 
is  driven  by  the  distribution  of  flow  angle  on  the  computational  boundary 
while  the  upstream  flow  is  driven  by  the  distribution  of  Riemann  variable 
associated  with  upstream  propagating  waves.  The  assumption  is  made  that  the 
far  field  flow  is  confined  by  parallel  walls,  which  is  not  restrictive  for 
duct  flow  problems. 


The  first-  and  second-  )rder  solutions  are  coupled  to  the  nonlinear 
numerical  solution  to  provide  a  smooth  transition  across  the  computational 
boundary  to  the  true  far  field  conditions  at  infinity.  The  coupling  is 
accomplished  by  the  boundary  conditions.  The  first-order  boundary  conditions 
provide  distributions  of  flow  quantities  to  be  imposed  along  the  boundary,  not 
constant  conditions.  The  second-order  boundary  conditions  provide  high 
accuracy  corrections  to  the  first-order  distributions.  The  first-  and 
second-order  boundary  conditions  represent  logical  asymptotic  extensions  of 
the  zero-order  (or  characteristic)  conditions.  Furthermore,  the  boundary 
analysis  can  be  coupled  with  any  inviscid  numerical  solution  method. 

Numerical  results  are  presented  for  a  duct/nozzle  configuration.  Results 
obtained  using  the  second-order  boundary  conditions  are  compared  with  those 
using  zero-  and  first-order  conditions.  An  extremely  fine  grid  was  used  for 
the  numerical  solution  in  order  to  accurately  assess  the  improvement  produced 
by  the  second-order  analysis.  It  was  found  that  the  computational  boundaries 
could  be  located  very  close  when  the  second-order  boundary  conditions  were 
used  with  no  loss  in  numerical  solution  accuracy.  The  reduced  size  of  the 
computational  field  further  reduced  the  number  of  grid  points  needed  for  the 
numerical  solution.  The  additional  computational  effort  required  to  upgrade 
the  boundary  conditions  to  second-order  is  very  slight.  The  benefits  of  the 
improved  boundary  conditions  are  also  achieved  when  a  coarser  grid  is  used. 


1 1 .   SECOND-ORDER  PERTURBATION  EQUATIONS 


The  system  of  two-dimensional,  steady,  linearized  Euler  equations  which 
describe  second-order  perturbations  from  a  constant  pressure  state  are  derived 
in  this  section.  For  reasons  outlined  in  Reference  1,  the  Riemann  variable 
formulation  of  Reference  2  will  again  be  used  because  of  its  close 
relationship  with  characteristic  (or  zero-order)  boundary  conditions  commonly 
used  in  numerical  solution  of  the  nonlinear  Euler  equations.  As  explained  in 
the  previous  section,  the  analysis  will  be  limited  to  isentropic  flow  in  ducts 
having  parallel  walls. 

The  two-dimensional,  isentropic  form  of  the  Euler  equations  is 
(Reference  2) 

^+  (q  -  a)  Vs  =   +  qa^  (2) 

89    _  8G     a2  8P  ,.,* 

It  +  q  5s  =  "  Vq  5H  (3) 

Velocity  magnitude  and  speed  of  sound  are  denoted  by  q  and  a,  respectively, 
and  P  is  the  logarithm  of  pressure.  The  Riemann  variables  Q  and  R  are  defined 
as 


Q  =  q  +  7TT  a 

R  =  q  -  — r  a 
Y-l 


(4) 


The  flow  angle  is  e,  time  is  denoted  by  t,  and  local  distances  along  and 
normal  to  the  streamline  direction  are  denoted  by  s  and  n,  respectively. 


For  steady  flow  the  analysis  can  be  greatly  simplified  by  defining  a  new 
dependent  variable 

T  =   Q  -  R  (5) 

and  replacing  equations  (1)  and  (2)  by 

<"2  -  »>  I  +  2  i  H  Is  ■  °  <6> 

a2  +  izl  q2  =  i  (7) 

The  local  Mach  number  is  denoted  by  M.  Equation  (6)  is  obtained  by 
subtracting  equations  (1)  and  (2).  Equation  (7)  is  obtained  by  adding 
equations  (1)  and  (2)  and  integrating.  The  constant  of  integration,  which  is 
proportional  to  stagnation  temperature,  can  be  set  to  unity  by  proper  choice 
of  non-dimensionalizing  quantities.  The  simplified  form  of  the  steady  Euler 
equations  is  then 

(M2  -  1)  |!+  2  q  H  §jj  =  0  (8) 

M2  99  +  _1_  I  31  =  o  (9) 

M  dS        y-1  T  3n   U  vy; 

a2  +  *±   q?  =  1  (10) 

In  regions  of  the  flowfield  where  nonlinear  effects  are  weak,  the  flow 
can  be  treated  as  a  perturbation  to  a  constant  pressure,  rectilinear  flow. 
Such  regions  occur  near  and  beyond  far  field  computational  boundaries.  The 
dependent  variables  in  equations  (8)  and  (9)  can  then  be  expanded  in 
asymptotic  series 

T  =  T^  +  Ti  +  T2  +  . . . 

(11) 

e  =  e„  +  ei  +  62  +  • • . 


The  flow  direction  at  infinity  is  assumed  constant  and  denoted  by  e^;  T^  is 
also  constant  because  the  flow  is  assumed  isentropic.  The  perturbation 
quantities  T-j  and  Q\   vanish  at  infinity.  For  duct  flow  between  parallel 
walls,  the  far  field  flow  angle  e^,  can  be  assumed  zero. 

Consistent  with  the  expansions  (11),  spatial  derivatives  in  equations  (8) 
and  (9)  can  be  approximated  by 


3    3  .  „  3    1  „2  3 

<  +  °2  8y-+  ••• 

(12) 


3S   8x  +  91  ^  "  2  Ql  3^  +  02  ay"  + 


3_  -  3_      3_   1  Q2  3_   Q  3_ 

3n  "  3y  "  8l  3x"  "  2  6l  3y  "  92  ax  +  ••• 

where  x  and  y  are  reference  Cartesian  coordinates  and  8i  and  02  are  measured 
from  the  x  axis  which  is  aligned  with  the  axis  of  the  duct  (see  Figure  1). 

If  expansions  (11)  and  (12)  are  introduced  into  equations  (8)  and  (9), 
the  resulting  first-order  perturbed  Euler  equations  are 

9  3T1        38i  ,  % 

B2  5T  -  2  «-M-  IT  '-   °  (13) 

39i    3Ti  ,   % 

The  second-order  equations  are 


9  3T2        8e2       a  r  3  2 
»2  3F  -  2  ^M~  W  -   ~   q^  [M~  91 


+  -4-  (1  +  *f±  \£   +  *f±  Mi)  Tf] 
2q„  Mro 


(15) 


2Q.M„^  +  ^  =  q^|ef-^-(l+^M!)T?]        (16) 


Velocity  and  Mach  number  at  infinity  are  constant  and  denoted  by  q,,,  and  Mm, 
respectively.  The  parameter  B  is  defined  by 


B  eVi-M^2  (17) 

Procedures  for  calculating  far  field  quantities  such  as  q*.  and  M,,,  are  outlined 
in  Reference  1. 


Asymptotic  expansions  of  the  Riemann  variables  Q  and  R  can  also  be 
defined  as 

Q  =  Q~  +  Qi  +  Q2  +  • • • 

(18) 
R  =  R^  +  Rt  +  R2  +  . .. 

Using  the  definition  (5)  and  the  expansion  (11)  for  T,  it  follows  that 

To,  =  Qc  -  R„ 

Ti  =  Qi  -  Ri  (19) 

T2  =  Q2  -  R2 


Introducing  the  expansions  (11)  and  (18)  into  the  algebraic  equation  (10) 
gives  the  first-  and  second-order  relationships 

(1+M„)  O4  -  (1-MJ  R!  =  0  (20) 

U+MJ  Q2  -  (1-MJ  R2  +     ""  .2  (1  +  *r   m£)  R?  =  0  (21) 

q<nl  l+"»J 

These  relations  will  be  used  later  in  Section  IV  where  the  boundary  conditions 
are  derived. 


Equations  (20)  and  (21)  can  also  be  used  to  express  the  perturbation 
equations  in  different  dependent  variables.  One  such  form,  which  will  be 
convenient  later  in  Sections  III  and  IV,  uses  R^  and  R2  in  place  of  T^  and  T2 
In  terms  of  these  variables,  the  first-order  equations  are 

v   3R1  891  ,     X 

(1-M.)  3/  +  q.  ^  =  0  (22) 

88i   3Ri 
q~  (1+M.)  9^  -  a/  =  0  (23) 

The  second-order  equations  are 

(1-M.)  |?2  +  Qa)  |£Z  =  1  Qm  JL  [M2  Q2  +  1    (1  +  Xxl  M^)  R?]   (24) 

3X      3y    '    3X         (1+M.)q£      * 

x  302   3R2   1    3  r  M~   2 

q.  (l+M.)  ar  -  ay  =  2  ^  ay  []Tm;  el 

(25) 
*— s  (1  +  2M.  +  *£   M3)  R2] 


(1+M„)2q£ 

The  second-order  equations  using  either  choice  of  dependent  variables  consist 
of  a  coupled  pair  of  non-homogeneous  equations  whose  right  hand  sides  (RHS) 
are  functions  of  the  first-order  solution  variables. 


III.   SO  .UTION  OF  SECOND-ORDER  EQUATIONS 


Solutions  of  the  second-order  perturbation  equations  for  duct  flows  are 
developed  in  this  section.  The  first-order  duct  flow  analysis  of  Reference  1 
showed  that  the  behavior  of  the  flow  in  the  downstream  far  field  is  governed 
by  the  distribution  of  flow  angle  e  on  the  downstream  computational  boundary. 
This  distribution  is  obtained  from  the  nonlinear  numerical  solution.  Also, 
the  first-order  downstream  boundary  conditions  for  the  computational  domain 
were  expressed  in  terms  of  the  perturbation  variable  Tj.  Therefore,  the 
second-order  equations  (15)  and  (16)  are  appropriate  for  analyzing  the  down- 
stream region. 

In  the  upstream  far  field  region,  the  flow  behavior  is  governed  by  the 
distribution  of  the  Riemann  variable  R  on  the  upstream  computational  boundary. 
This  distribution  is  obtained  from  the  nonlinear  numerical  solution.  The 
first-order  upstream  boundary  conditions  for  the  computational  domain  were 
expressed  in  terms  of  the  perturbation  variables  ei  and  Qj.  Therefore,  the 
second-order  perturbation  equations  (21),  (24)  and  (25)  are  appropriate  for 
analyzing  the  upstream  region. 

Downstream  Region 

The  second-order  equations  (15)  and  (16)  can  be  combined  into  the  single 
equation 


2  !f!2  +  ^  .  JL   [M£e?  +  *  (1  +  iji  mZ)  if]  (26) 

8x2    9y2    8x3y        4q£      2 


From  Reference  1,  the  solution  of  the  first-order  system  (13)  and  (14)  for 
duct  flows  can  be  expressed  as 


6!   =  z  An  sin  nny  e"nil)(/B 
1 

(27) 

^qJL  -  -mrx/B 

Tj     = —  E  An  cos  riTiy  e 

*       1 


where  An  represents  the  Fourier  coefficients  of  the  distribution  of  flow  angle 
e  on  the  downstream  computational  boundary  located  at  x=0.  The  width  of  the 
duct  has  been  taken  as  unity.  Note  that  the  solution  (27)  for  Ti  does  not 
include  a  zero-mode  term  (i.e.,  Aq)  which  is  related  to  the  mean  value  of  the 
perturbation  T-T^  at  the  boundary.  It  was  pointed  out  in  Reference  1  that 
such  a  term  is  a  second-order  effect.  This  will  be  verified  below.  Using  the 
first-order  solution  in  the  RHS,  equation  (26)  becomes 


~     829?  329?  TT2M0D    co    a>  o         p 

Bz  — ^r   +  — ~r   =  5  l  l   Am  An  W   +  "V  M~)  (mz-nz)  sin(m-n)ny 

Bx^    9y^    26J  i  i  * 

(28) 
+  ^  \l   (m+n)2  sin(m+n)T,y]  e-(m+n)7TX/B 


To  satisfy  the  duct  wall  boundary  conditions  of  zero  flow  angle,  a 
particular  solution  of  equation  (28)  must  have  the  form 


82  =  t~T  *  *■   Hmn(*)  sin(m±n)Tiy  (29) 

26J  i  i 


This  results  in  the  two  component  ordinary  differential  equations  (depending 
on  choice  of  sign)  for  each  value  of  m  and  n 

62  Hmn  -  ,2(m+n)2  Hmn  =  X+l  Mz  (m+n)2  A(n  Ar  e-(m+n)nx/B  (3Q) 

BZ  H^n  -  n2(m-n)2  Hmn  ■  (2  +  *£   m£)  (l2-n2)  Am  An  e^m+n^x/B     (31) 
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The  respective  solutions  of  these  equations  are 
.2 


Mnn 


=  -  *i>  £  (m+n)  Am  A„  x  t-l™)«l*  (32) 


mn  '  77  (    2   <"'  Urn        m     n  (33' 


4tt< 


Referring  to  the  first-order  solution  (27),  the  solution  for  ei  can  be 
used  to  exactly  represent  the  computational  boundary  distribution  of  flow 
angle  obtained  from  the  nonlinear  numerical  solution.  This  requires  that  92 
vanish  at  x=0.  This  can  be  enforced  by  adding  a  solution  of  the  homogeneous 
portion  of  equation  (31)  to  the  solution  (33).  The  full  second-order  solution 
for  e  in  the  downstream  region  is  then 


e  =  l   An  sin  mry  e"nitx/B  +  k{   x  1  1   (m+n)  Am  An  sin(m+n)*y  e"(m+n)l,x/B 


1  1  1 


(34) 


+  k2  J  J  ^  A,,  A„  sin(m-n)*y  frl™)**'*   -  e- 1 m-n | ^x/B -, 

1  1  mn 


where 


kl  E  -  8  04 


Y+l  ^M0 


(35) 


The  solution  for  T2  can  be  obtained  from  equation  (15)  using  the  solution 
(29)  for  82  which  was  determined  above.  After  integrating,  the  full 
second-order  solution  for  T  in  the  downstream  region  is 
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t       t         2qmM°°  £  n     „„P   „  „  «-nnx/B 
T  =  Tm  -  — - —  E  An  cos  nny  e 

14       1 


+  Ki  E  E  Am  An  cos(m+n)iiy  [1  +  J  (m+n)x]  e-(m+n)7TX/B 
1   1  B 

(36) 

K2n^  cos(m-n),y  [(m-n)2  e-(<"+n),x/B  _    |fn2_n2|  e-\m-n\.x/^ 


+  E  E  Am  An  [K3  cos(m-n)Tiy  +  K4  cos(m+n)Tiy]  e-(m+n)ltX/B 
1  1 


The  coefficients  are  defined  by 


K    1+1  ^ 
Kl  =  4   «4 


5 
i 


s  1  "-4   (2  +  ^  Hi) 


K2  "4  34 


(37) 


qJJ»  .         y+1  m2  ,  Y-3  m4\ 
<3  =  "  -74-  (1  +  4  M»  +  4  M„) 


Y  ^mM°°  M,  I"3  M2     l±i  M41 

K4  =  -  -^-  (1+  ^4-  M.  +  -^-  Mo,) 


The  solution  (36)  contains  longitudinally  decaying  plane  waves 
corresponding  to  m=n.  They  are  represented  by 


T2  =  K3  E  An  e 
1 

At  x=0  the  quantity  T2  represents  the  mean  value  of  the  perturbation  T-Tm  at 
the  boundary. 


(38) 
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Upstream  Region 

The  second-order  equations  (24)  and  (25)  can  be  combined  into  the  single 
equation 

9  3^R9   32R?   i     a2    ?        ?   i       v  1   "*   ? 

6  ^r  +  rr    2  q~  (!"2  [M~  (1+M~>  ei +  S  <*  +  ?  M~>  Ri] 

3x^    3y^    '     8x^  ,/       t 

(39) 

-  &  ^ -l "  (wfe?  <»  *  »*■  ♦  ^  *  *' 


From  Reference  1,  the  solution  of  the  first-order  system  (22)  and  (23)  can  be 
expressed  as 

R!  =  z   Cn  cos  nny  enilX/B 
1 


(40) 


ei =  uiky  J c" sin  nny  e"nx/B 


where  Cn  represents  the  Fourier  coefficients  of  the  distribution  of  the 
Riemann  variable  perturbation  R-R^  on  the  upstream  computational  boundary 
located  at  x=0.  As  explained  in  Reference  1,  the  mean  value  of  R-R^  on  the 
boundary  represented  by  the  zero-mode  coefficient  Cq  is  a  second-order  effect 
and  is  not  included  in  the  first-order  solution  (40).  Using  the  first-order 
solution  in  the  RHS,  equation  (39)  becomes 

2  32R2   a2R2  ..  »?    l    rn  +  m  +  m2  +  id  m3  +  id  mS 

6     372"  +  3^2"  "   4^2     q.(l+M.)   C{1   +  M"  +  M"  +     2     M~  +     4     M-> 

n  Cm  Cn  (m+n)2  cos(m-n)uy 

1    1  (41) 

-   (1  +  M*,  -  3m£  +  *d  M3  _  id  Mf)  J  J  Cfn  Cn   (m-n)2  cos(m-n)Tiy 


1   1 


+  (y+1)  m!  l  E  Cm  Cn  (m+n)2  cos(m+n)iiy]  e(m+n),,x/6 
1   1 
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Following  the  solution  procedure  used  for  the  downstream  region,  a 
particular  solution  of  equation  (41)  must  have  the  form 

R2  =  7     "     E  l   Gmn(x)  cos(m±n)Try  (42) 


This  results  in  the  two  component  ordinary  differential  equations  (depending 
on  choice  of  sign)  for  each  value  of  m  and  n 

62  Gmn  -  *2(m+n)2  G(mi  .  {r+l)  £   (m+n)2  Cfn  Cn  e(m+n)lTX/e  (43) 

B2  G^n  -  -2(m-n)2  Gmn  =  [b^m+n^  +  b2(m-n)2]  cm  Cn  e(m+n)l,x/B     (44) 
The  coefficients  b^  and  b2  are  defined  by 

bi  =  1  +  Mm  +  m£  +  ^  M^  +  *jr  mJ 

b2  =  -  (1  +  M.  -  3m£  +  ^  M^  -  ^  m!) 


(45) 


The  respective  solutions  of  these  equations  are 

4 
Gm„  -  ^  £  («")  Cm  Cn  x  e<-n>"/'J  (46) 

6mn  ■  ^2  t">l  I™-)2  -  b2  (m-n)2]  ^  ^  e<m+n>"</8  (47) 


Referring  to  the  first-order  solution  (40),  with  the  exception  of  the 
mean  value,  the  solution  for  R^  can  be  used  to  exactly  represent  the 
computational  boundary  distribution  of  R-Rm  obtained  from  the  nonlinear 
numerical  solution.  This  requires  that  all  Fourier  modes  of  R2  except  the 
zero  mode  vanish  at  x=0.  This  can  be  enforced  by  adding  a  solution  of  the 
homogeneous  portion  of  equation  (44).  The  solution  for  R2  in  the  upstream 
region  is  then 
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R2  =  K5  x  z  I  (m+n)  Cm  Cn  cos(m+n)„y  e(m+n)7,x/B  +  4K6  i  C2n  e2nnx/B 
1   1  1 

H6H  ^  Cm  Cn  cos(m-n)*y  [e*^*/*  -  e|m-n|.x/63 

1   1       mn 


(48) 


H7H  ^^  Cm  Cn  cos(m-n)»y  [e^")"^  -  elm-ni"x/B] 
1   1       mn 


where 

K5  E 


BB3q.(l+H.) 

K6  =  7t4—  (1   +  m£  +  *z3  m2)  (49) 

16B*q.  * 


The  corresponding  solution  for  Q2  can  be  obtained  from  the  second-order 
relation  (21)  as 

l"Ma>  Mc  v_1         9  9 

02  =  rat  "z "  qJT^P  <»  ♦  V  H->  «f  <5°> 


The  solution  for  e2  can  be  obtained  from  equation  (24)  using  the  solution 
(48)  for  R2  which  was  determined  above.  After  integrating,  the  full 
second-order  solution  for  e  in  the  upstream  region  is 


: 
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9  =  n  m?m  a  ^  cn  sin  rmy  en"x/B  +  k3  I  z  Cm  Cn  sin(m+n)iiy  e(m+n)t,x/B 

Moo  i.  1+rlco  J       1  11 

+  k4  x  E  z  (m+n)  Cm  Cn  sin(m+n)*y  e(m+n)nx/B 

(51) 

L     r  r  m+n  r     r     .jH/n  „\       r/m  -\  |m-n |      lm-n|iix/0       /m     *   a(m+n)iix/Bi 
+  k5  E  E  "in  cm  cn  sin(m-n)*y  t(m+n)-Ljjj^J-  e1       '     /p  -  (m-n)  ev  J 

m/n 

H6n^cmcn  sin(m-n)ny  [(m+n)  e<m+n>"x/e  -  |m-n|  el"""!""] 
l  l  mn 


The  coefficients  are  defined  by 


k3  =  ^ o  [B2  (1  +  Mro)  +  *+I  m^ 

4eq£(l+Mro)2 

.4 


k4  = 


(+1  TiMo 


8     „:>  2 


B2q£(l+Ma>)2 


(52) 


k5  e  ^ (1  +  m£  +  ^  m£) 

16Bq£(l+Mm) 


16Bq£(l+Mro)2 


(1  +  Mm  -  3m£  +  *|I  mI   -  *£■   Ml) 


The  solution  (48)  also  contains  longitudinally  decaying  plane  waves 
corresponding  to  m=n.  They  are  represented  by 


1 


1   (1  +  „£  +  t2  ft  -  c2  e2n«/B 

46^0.  * 


(53) 

r2  2nwx/B 
Ln  e 

1 
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From  equation  (50)  the  corresponding  plane  waves  for  Q  are 


^  =  A     r?  M  13  (1  '  M~  +  M~  "  *T  M"  +  ¥  M~>  =  C"  e2n7,X/B      (54) 

"Mco  V  A  '"''00  j 


At  x=0  the  quantities  Q?  and  R?  represent  the  mean  values  of  the  perturbations 

Q-Qro  and  R-Rro  at  the  boundary,  respectively.  Note  that  the  mean  value  of  the 

R  perturbation  is  predicted  by  the  analysis  and  not  determined  directly  from 
the  numerical  solution  as  are  the  other  Fourier  modes. 
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IV.   BOUNDARY  CONDITIONS  DEVELOPMENT 


Second-order  boundary  conditions  for  duct  flows  are  developed  in  this 
section  based  on  the  linearized  Euler  solutions  obtained  in  the  previous 
section.  They  represent  a  logical  extension  of  the  first-order  boundary 
conditions  developed  in  Reference  1. 

Downstream  Boundary 

The  downstream  boundary  conditions  are  derived  from  the  definition  (5) 
rewritten  as 

R  =  Q  -  T  (55) 

Using  the  solution  (36)  for  T  and  applying  the  relation  (55)  at  the  boundary 
(assumed  located  at  x=0),  the  distribution  of  R  on  the  boundary  (i.e.,  the 
boundary  conditions)  is  calculated  according  to 


4      ZqooMco  °° 
Rb  =  Qnum  -  ^ry  a~  +  ~j~    *   An  cos  nr,y 


+  (K\   +  K4)  EE  Am  An  cos(m+n)ny  (56) 

1  1 


+  ;  ;  [K  C-n)2-  |m2-n2|  +  ^  ^  An  cos(m_n)„. 


1  1  mn 

The  boundary  distribution  of  Q  obtained  from  the  nonlinear  numerical  solution 
is  denoted  by  Qnum-  TnG  zero  Courier  mode  corresponding  to  m=n  in-  this 
expression  is  proportional  to  the  coefficient  K3.  It  represents  the  mean 
value  of  the  R  perturbation  at  the  boundary  and  is  clearly  a  second-order 
effect. 
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Upstream  Boundary 

The  upstream  boundary  conditions  are  derived  by  applying  the  perturbation 
relation  (50)  and  the  solutions  (48)  and  (51)  at  the  boundary  (assumed  located 
at  x=0).  From  these  relations  the  distributions  of  Q  and  e  on  the  boundary 
(i.e.,  the  boundary  conditions)  are  calculated  according  to 

1    _M  CO  CO  O 

Qb  =  Q~  +  T7IT  f^  cn  cos  n7,v  +  4  K6  E  Cn] 

1    •  ™CO  1  1 


K 


-=  (1  +  *j±  m£)   [e  Cn  cos  n*y]2 
qco(l+Mco)J  c  i 


(57) 


and 


eb  =       m+mj     E  cn  sin  ni,y  +  k3  l  E  cm  cn  sin(m+n)ny 

+  k5  "     ST  l1*^  -<m-n>]  C*  C"  sin(m-n)ny  (58) 

ntfn 


+  k6  E  E  JHr  [m+n  -  |m-n|]  Cm  Cn  sin(m-n)ny 
1  1  mn 

The  zero  Fourier  mode  proportional  to  K6  in  equation  (57)  is  again  evident. 
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V.   APPLICATIONS 


Numerical  solutions  of  the  Euler  equations  have  been  calculated  for 
two-dimensional,  isentropic,  steady  duct/nozzle  flow  using  the  second-order 
boundary  condition  procedures  developed  in  the  previous  section.  Second-order 
boundary  condition  results  are  compared  with  those  produced  using  the  conven- 
tional zero-order  characteristic  boundary  conditions  and  also  the  first-order 
boundary  conditions  developed  in  Reference  1. 

A  solution  algorithm  was  used  for  the  nonlinear  Euler  equations  (1)  -  (3) 
which  is  based  on  the  method  presented  in  Reference  2.  It  uses  explicit  time 
integration  to  relax  to  steady  state  conditions.  It  should  be  emphasized  that 
the  boundary  condition  analysis  is  independent  of  the  choice  of  inviscid, 
nonlinear  solution  method. 

The  duct/nozzle  geometry  is  shown  schematically  in  Figure  1.  The  flow  is 
characterized  by  pro,  the  downstream  pressure  at  infinity,  which  produces  a 
mass  flow  per  unit  area  w  through  the  duct.  The  linearized  solutions  given  by 
equations  (34),  (36),  (48),  (50)  and  (51)  are  assumed  valid  in  the  semi- 
infinite  regions  I  and  III  and  the  computational  boundary  conditions  are 
applied  at  the  upstream  and  downstream  boundaries  AA  and  BB  of  the  nonlinear 
computational  region  II. 

The  actual  shape  of  the  duct/nozzle  and  the  computational  grid  are  shown 
in  Figure  2.  The  shape  is  identical  to  that  used  for  the  calculations  pre- 
sented in  Reference  1.  The  nozzle  contour  is  sinusoidal  and  symmetric  about 
the  centerline.  The  computational  grid  for  this  portion  of  the  nozzle  had 
dimensions  81  x  41,  which  is  twice  the  size  of  that  used  in  Reference  1.  The 
increased  grid  dimensions  allowed  second-order  effects  to  be  quantified  more 
accurately.  The  area  ratio  of  the  nozzle  is  .75  and  the  upstream  and  down- 
stream areas  are  equal.  For  these  constant  area  sections  of  the  duct,  addi- 
tional rectangular  grid  cells  could  be  added  without  altering  the  basic  81  x 
41  grid.  This  served  to  minimize  the  effect  of  grid  changes  on  the  calcula- 
tions when  the  computational  boundaries  were  moved  in  order  to  assess  the 
accuracy  of  the  boundary  conditions. 
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Because  the  configuration  is  symmetric,  calculations  were  limited  to  the 
lower  half  of  the  duct  and  .1  centerline  symmetry  condition  was  used.  Although 
the  configuration  used  for  these  calculations  is  simple,  the  boundary  condi- 
tion analysis  of  the  previous  section  is  general  and  applicable  to  unsymmetric 
configurations  having  unequal  upstream  and  downstream  areas.  Use  of  the 
simple  configuration  is  sufficient  to  demonstrate  the  validity  of  the 
analysis. 

The  second-order  upstream  and  downstream  boundary  conditions  are  given  by 
equations  (56),  (57)  and  (58).  The  associated  analytic  far  field  solutions 
are  given  by  equations  (21),  (34)  and  (36)  for  the  downstream  region  and  by 
equations  (48),  (50)  and  (51)  for  the  upsteam  region.  The  zero-order  (or 
characteristic)  boundary  conditions  consist  of  imposing  the  constant  value  of 
()„  and  a  zero  value  of  e  along  the  upstream  boundary  and  the  constant  value  of 
Rm  along  the  downstream  boundary.  The  first-order  upstream  boundary  condi- 
tions consist  of  imposing  additional  distributions  of  Qi  and  ei  along  the 
boundary  as  determined  from  equations  (20)  and  (40).  The  first-order  down- 
stream boundary  conditions  consist  of  imposing  an  additional  distribution  of 
Rl  along  the  boundary  as  determined  from  equations  (27)  and  (55).  Additional 
details  are  given  in  Reference  1. 

Computational  results  are  presented  for  a  single  value  of  Po,  but  with 
the  computational  boundaries  located  at  two  different  longitudinal  stations. 
One  boundary  is  sufficiently  far  removed  so  that  all  three  sets  of  boundary 
conditions  produced  essentially  the  same  results  within  the  computational 
domain.  The  second  boundary  is  extremely  close  to  the  nozzle  portion  of  the 
duct  so  that  the  relative  accuracy  of  the  various  boundary  conditions  can  then 
be  evaluated. 

Results  obtained  using  the  complete  grid  shown  in  Figure  2  are  presented 
in  Figures  3  and  4.  This  grid  has  40  columns  of  grid  cells  in  both  the 
upstream  and  downstream  constant  area  portions  of  the  duct.  As  mentioned 
above  the  results  were  nearly  identical  for  all  three  sets  of  boundary  condi- 
tions. Figure  3  shows  pressure  and  Mach  number  distributions  along  the 
centerline  and  lower  wall  of  the  duct/nozzle.  Pressure,  Mach  number,  and  flow 
angle  contours  are  presented  in  Figure  4.  These  results  serve  as  a  reference 
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for  evaluating  the  accuracy  of  solutions  where  the  computational  boundaries 
are  moved  closer  to  the  nozzle  portion  of  the  duct. 

Results  for  the  shortened  computational  domain  obtained  using  the  zero- 
order  boundary  conditions  are  presented  in  Figures  5  and  6.  There  was  only 
one  column  of  grid  cells  in  both  the  upstream  and  downstream  constant  area 
portions  of  the  duct  for  this  case.  Pressure  and  Mach  number  distributions 
are  shown  in  Figure  5  and  contours  are  shown  in  Figure  6.  The  degradation  in 
the  solution  as  a  result  of  the  boundary  proximity  is  very  evident. 

Application  of  the  first-order  boundary  conditions  to  the  short  computa- 
tional domain  gives  the  results  shown  in  Figures  7  and  8.  Eleven  Fourier 
modes  were  needed  to  accurately  describe  the  boundary  distributions  of  e  and  R 
for  this  extremely  close  boundary  location.  Pressure  and  Mach  number 
distributions  are  shown  in  Figure  7  and  contours  in  Figure  8.  Linearized 
first-order  solution  results  obtained  from  equations  (20),  (27)  and  (40)  have 
been  added  upstream  and  downstream  of  the  computational  boundaries. 
Significant  improvement  is  provided  by  the  first-order  boundary  conditions. 
Some  mismatch  at  the  boundary  is  evident,  but  this  can  be  attributed  to 
neglect  of  second-order  interactions. 

Finally,  results  for  the  shortened  computational  domain  using  the 
second-order  boundary  conditions  are  presented  in  Figures  9  and  10.  Eleven 
Fourier  modes  were  also  used  in  obtaining  these  results.  Pressure  and  Mach 
number  distributions  are  shown  in  Figure  9  and  contours  in  Figure  10.  The 
results  within  the  numerical  solution  portion  of  the  domain  are  nearly 
identical  to  those  shown  in  Figures  3  and  4.  Linearized  second-order  solution 
results  obtained  from  equations  (48),  (50)  and  (51)  and  from  equations  (21), 
(34)  and  (36)  have  been  added  upstream  and  downstream  of  the  computational 
boundaries,  respectively.   It  is  evident  that  the  linearized  far  field 
analytic  solutions  provide  for  a  smooth  transition  across  the  computational 
boundary  to  the  true  far  field  conditions  at  infinity. 

A  more  quantitative  comparison  between  the  three  sets  of  boundary  condi- 
tions can  be  obtained  by  examining  the  distribution  of  flow  variables  along 
the  computational  boundaries.  Figure  11  shows  a  comparison  of  the  distribu- 


22 


tion  of  pressure  and  flow  .  ngle  along  the  upstream  and  downstream  boundaries 
for  the  short  computational  domain.  Results  at  the  same  longitudinal  loca- 
tions taken  from  the  reference  numerical  solution  using  the  long  grid  (Figure 
3)  are  also  shown.  These  results  exhibit  classical  asymptotic  convergence  to 
the  reference  solution.  The  larger  increment  between  the  zero-order  and 
first-order  pressure  results  at  the  upstream  boundary  is  primarily  due  to  the 
fact  that  the  zero-order  boundary  conditions  impose  a  zero  flow  angle  along 
the  boundary. 

It  was  pointed  out  in  Section  III  and  also  in  Reference  1  that  the  mean 
values  of  the  perturbations  at  the  computational  boundaries  are  neglected 
within  the  first-order  analysis.  In  Section  III  they  were  shown  to  appear  as 
second-order  effects  described  by  equations  (38)  and  (53)  evaluated  at  x^O. 
In  the  downstream  region  equation  (38)  indicates  that  the  relationship 


00  2-,-1  f1 
■  E  An3 

1  J  o 


(T-T.)  dy»  -5J=(1  *3t+i  m£  +  ^h!)  (59) 


should  hold.  Likewise,  from  equation  (53)  the  corresponding  relation  for  the 
upsteam  region  is 


[E  CnY1  f  (R-Rco)  dy  =  -±—   (1  +  m£  +  *£■  m2) 

l         J  o  4^q«  d 


(60) 


These  relations  are  shown  in  Figure  12  for  P™  =  0.90.  The  Fourier 
coefficients  and  integrals  are  determined  at  various  longitudinal  upstream  and 
downstream  locations  using  flowfield  data  taken  from  the  reference  long  grid 
solution  of  Figure  3.  Distance  measured  from  the  start  of  the  constant  area 
portions  of  the  duct  is  denoted  by  £.  These  results  confirm  the  second-order 
nature  of  the  perturbation  mean  values. 
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It  can  be  seen  from  equations  (56),  (57)  and  (58)  that  very  little 
additional  computational  effort  is  required  to  upgrade  the  boundary  conditions 
from  first-order  to  second-order.  Most  of  the  computational  effort  at  the 
boundaries  is  spent  determining  the  Fourier  coefficients  which  are  needed  to 
impose  first-order  boundary  conditions.  Therefore,  computational  boundaries 
having  isentropic  inflow  can  be  upgraded  easily  even  though  a  source  of 
entropy  production  may  exist  within  the  nonlinear  computational  domain. 

The  results  presented  in  Figures  3  through  12  were  obtained  using  a  very 
fine  computational  grid  in  order  to  accurately  quantify  the  differences 
between  the  boundary  condition  models.  To  demonstrate  that  the  improvement 
provided  by  the  second-order  procedure  can  also  be  achieved  with  a  moderately 
sized  grid,  the  calculations  were  repeated  with  the  grid  used  in  Reference  1. 
One  column  of  grid  cells  was  used  in  both  the  upstream  and  downstream  constant 
area  portions  of  the  duct.  These  results  are  shown  in  Figures  13  and  14. 
Pressure  and  Mach  number  distributions  are  shown  in  Figure  13.  Distributions 
of  pressure  and  flow  angle  along  the  upstream  and  downstream  computational 
boundaries  are  compared  in  Figure  14.  Note  that  these  boundaries  are  located 
slightly  upstream  and  downstream,  respectively,  from  those  of  the  fine  grid 
solution.  Asymptotic  convergence  to  the  reference  solution  (Reference  1)  is 
again  evident.  Eleven  Fourier  modes  were  needed  to  resolve  the  boundary 
distributions  of  e  and  R. 
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VI.   SUMMARY 


Highly  accurate  second-order  far  field  boundary  conditions  for 
two-dimensional,  isentropic  duct  flow  have  been  developed  from  analytic 
solutions  of  the  second-order  linearized  Euler  equations.  These  boundary 
conditions  represent  a  logical  extension,  in  an  asymptotic  sense,  of  the 
first-order  boundary  conditions  derived  in  Reference  1  from  analytic  solutions 
of  the  first-order  linearized  Euler  equations.  The  first-order  boundary 
conditions  of  Reference  1  are  themselves  a  logical  extension  of  the  zero-order 
(or  characteristic)  boundary  conditions  commonly  used  in  numerical  solutions 
of  the  nonlinear  fluid  dynamic  equations.  The  boundary  conditions  and 
analytic  solutions  provide  a  smooth  transition  across  a  computational  boundary 
to  the  true  far  field  conditions  at  infinity.  The  boundary  procedure  is 
general  in  that  it  can  be  used  in  conjunction  with  any  numerical  solution 
method. 

The  second-order  linearized  Euler  equations  are  non-homogeneous  and  have 
been  solved  exactly  using  separation  of  variables  and  Eourier  analysis.  This 
procedure  was  used  in  Reference  1  to  solve  the  linearized,  homogeneous 
first-order  equations.  Extension  of  the  second-order  analysis  to  allow 
non-isentropic  flow  conditions  is  probably  not  justified  because  solution  of 
the  non-isentropic  equations  requires  an  approximate  two-level  perturbation 
procedure  (see  Reference  1  for  details).  However,  even  if  the  flow  being 
analyzed  has  an  entropy  source  in  the  nonlinear  region,  the  second-order 
boundary  conditions  can  still  be  applied  at  the  inflow  computational  boundary 
if  the  incoming  flow  is  isentropic. 

Use  of  zero-order  (or  characteristic)  boundary  conditions  requires  that 
the  computational  boundaries  be  located  far  from  the  nonlinear  region  of  the 
flow.  Closer  placement  of  the  boundaries  may  result  in  a  significant  amount 
of  solution  degradation.  The  first-order  boundary  conditions  derived  in 
Reference  1  allow  the  boundaries  to  be  located  much  closer  thereby  reducing 
the  number  of  grid  points  needed  for  the  numerical  solution  and  also  the 
number  of  iterations  for  solution  convergence.  This  leads  to  a  significant 
reduction  in  the  amount  of  computational  effort  required  for  the  nonlinear 


2b 


numerical  solution  because  the  additional  calculations  required  for  the 
first-order  boundary  conditions  is  modest. 

Use  of  the  second-order  boundary  conditions  allows  the  computational 
boundary  to  be  placed  even  closer  with  a  further  reduction  in  the  number  of 
grid  points.  The  amount  of  additional  computational  and  implementation  effort 
is  very  slight  so  that  isentropic  inflow  boundary  conditions  can  be  upgraded 
with  little  penalty. 

In  Reference  1  it  was  pointed  out  that  the  mean  values  of  the  flow 
variable  perturbations  at  the  computational  boundary  were  not  described  within 
the  first-order  analysis.  It  was  further  postulated  that  these  mean  values 
represented  a  second-order  effect.  The  nature  of  this  second-order  interac- 
tion was  clarified  by  the  present  analysis  and  verified  by  numerical  results. 
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Figure  1.    Duct/Nozzle  Schematic 


Figure  2.  Duct/Nozzle  Geometry  and  Grid 
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Figure  3.    Pressure  and  Mach  Number  Distributions 
Zero-,  First-  and  Second-Order  Boundary  Conditions 
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Figure  4.    Pressure,  Mach  Number  and  Flow  Angle  Contours 
Zero-,  First-  and  Second-Order  Boundary  Conditions 
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Figure  5.    Pressure  and  Mach  Number  Distributions 
Zero-Order  Boundary  Conditions 
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Figure  6.    Pressure,  Mach  Number  and  Flow  Angle  Contours 

Zero-Order  Boundary  Conditions 
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Figure  7.    Pressure  and  Mach  Number  Distributions 
First-Order  Boundary  Conditions 
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Figure  8.    Pressure,  Mach  Number  and  Flow  Angle  Contours 

First-Order  Boundary  Conditions 
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Figure  9.    Pressure  and  Mach  Number  Distributions 
Second-Order  Boundary  Conditions 
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Figure  10.    Pressure,  Mach  Number  and  Flow  Angle  Contours 
Second-Order  Boundary  Conditions 
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Figure  11.  Pressure  and  Flow  Angle  Distributions 

on  Computational  Boundaries 
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Figure  13.    Pressure  and  Mach  Number  Distributions 
Second-Order  Boundary  Conditions 
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Figure  14.    Pressure  and  Flow  Angle  Distributions 

on  Computational  Boundaries 
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